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Abstract 

It is shown that tunnelling splittings in ergodic double wells and resonant widths in ergodic 
metastable wells can be approximated as easily-calculated matrix elements involving the wavefunc- 
tion in the neighbourhood of a certain real orbit. This orbit is a continuation of the complex orbit 
which crosses the barrier with minimum imaginary action. The matrix element is computed by 
integrating across the orbit in a surface of section representation, and uses only the wavefunction 
in the allowed region and the stability properties of the orbit. When the real orbit is periodic, the 
matrix element is a natural measure of the degree of scarring of the wavefunction. This scarring 
measure is canonically invariant and independent of the choice of surface of section, within semiclas- 
sical error. The result can alternatively be interpretated as the autocorrelation function of the state 
with respect to a transfer operator which quantises a certain complex surface of section mapping. 
The formula provides an efficient numerical method to compute tunnelling rates while avoiding the 
need for the exceedingly precise diagonalisation endemic to numerical tunnelling calculations. 


1 Introduction 

Tunnelling in one dimension is an extremely well-studied and understood problem jij. Using standard 
WKB methods, together with appropriate uniform approximations, one can accurately approximate 
such quantities as subbarrier transmission amplitudes, energy level splittings in double wells and reso¬ 
nance lifetimes in quasi-bound systems. In general this is not true in higher dimensions. The problem 
is often enormously more difficult, for the simple reason that the underlying classical mechanics is 
qualitatively richer. 

This situation is especially evident in problems with classically chaotic motion. In such cases, there 
is no simple constructive theory for wavefunctions and any incorporation of tunnelling effects must 
confront this fact. It is certainly impossible to write simple closed expressions for a tunnelling rate 
associated with an individual state along the lines of what exists in one dimension. In this paper we 
develop a formula which comes as close as possible within these limitations. We concentrate on energy 
level splittings between states in a symmetric double well, though the calculation for resonance widths 
is also summarised. Our calculation is built upon an approach developed by Auerbach and Kivelson 
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||] and by Wilkinson and Hannay ||], henceforth referred to as AK and WH respectively. They use 
semiclassical Green’s functions, constructed from complex trajectories, to extend wavefunctions from 
allowed to forbidden regions and thence to calculate spectral tunnelling effects. We show how this 
method can be used to express the splitting as a matrix element which measures the concentration of 
the wavefunction in the neighbourhood of a particular classical trajectory emerging from the optimal 
tunnelling route. 

The formula assumes that the wavefunction is known (possibly following numerical diagonalisa- 
tion), but otherwise uses quantities that are straightforwardly calculated from classical dynamics. 
The technical benefit is that, while a completely numerical determination of the splitting requires 
that diagonalisation be performed to a precision at least comparable with the energy level splitting 
(which can be extreme), this matrix element requires only that the wavefunction be calculated to a 
precision comparable with standard semiclassical errors. Perhaps more importantly, there is a theo¬ 
retical benefit in that tunnelling rates are simply and directly related to the the morphology of the 
wavefunction in the allowed region. This can be used as a starting point for statistical analyses and 
also for semiclassical analysis using Green’s functions etc. 

When the trajectory underlying the matrix element is periodic, the splitting is nothing other than 
a scarring weight of the wavefunction. Scarring Q, or the excessive accumulation of certain eigen¬ 
states around periodic orbits, has become a common means of characterising chaotic wavefunctions. 
Intuitively, one expects exceptionally high tunnelling rates to occur in states having a large ampli¬ 
tude along the optimal tunnelling route. The matrix element quantifies this idea. It measures the 
amplitude of the wavefunction along a certain real trajectory which connects to what we will call 
the bounce orbit — an imaginary-time complex trajectory which crosses the barrier with minimal 
imaginary action. One often finds that the bounce orbit lies on a symmetry axis, in which case its 
real extension is periodic and the terminology of scarring is appropriate. The formula works even if 
the orbit is not periodic, however, though the scarring interpretation is then weaker. 

There has recently been a lot of interest in the current-voltage characteristics of tunnelling diodes 
in the presence of crossed electric and magnetic fields H Hi- I n those systems it has been argued 
that the tunnelling current is strongly mediated by scarred states. In this paper we quantify that 
very connection between scarring and tunnelling but in contexts where the tunnelling segment plays 
an active role in dynamics (up to now, tunnelling diodes have been treated assuming the barriers 
to be sufficiently thin that tunnelling can be incorporated simply using reflection and transmission 
coefficients). There are some formal similarities between our results and certain calculations for 
tunnelling diodes — our configuration space formulation to that of |7j and the phase space formulation 
to that of Q. There are also differences however. We assume the wave function is highly excited and 
fluctuating strongly and get a matrix element with a strongly localised Green’s function. In |7], |8| 
one finds the contrary: the wave function is localised because incoming electrons are assumed to be 
in an approximate ground state and the Green’s function, associated with real dynamics, is strongly 
oscillatory. 

The paper is divided as follows. An important element in getting simple and universal expressions 
for the splitting is to normalise it by an average background associated with the bounce orbit. In Q 
we derived an expression for the average splitting. Controlled by the imaginary action of the bounce 
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orbit, this average varies strongly with energy and with system parameters. Dividing each splitting by 
the average splitting at that energy leaves a set of rescaled splittings which can usefully be compared 
between states lying in very different parts of the spectrum. We begin in section || by reviewing the 
most important elements of this calculation. 

In section || we review the calculations of AK and WH and give a derivation of the main matrix 
element. The review is necessary because certain technical aspects of the calculation, having to do 
with choice of Green’s functions etc., must be stressed in order to understand our result. In addition 
we implement the formula somewhat differently; we consider the value of the wavefunction along a 
cut inside the well rather than on the energy contour as done in WH or along a cut in the forbidden 
region as in AK. Having restated the result, an interpretation is developed as a matrix element in a 
Poincare-section representation similar to that used in the transfer operator method of Bogomolny 
[10]. In a Wigner-Weyl calculus restricted to the Poincare section, the rescaled splitting is simply 
the overlap of the Wigner function of the state with a Gaussian centered on the real extension of the 
bounce orbit. The exponent of this Gaussian is directly related to the complex stability matrix along 
this tunnelling orbit and is a quadratic function on the Poincare section. In section |5] we outline how 
these results may be translated to the calculation of resonance widths in metastable wells. In section ||| 
the general results are given a formal interpretation in terms of transfer operators and this is then 
used as the basis for a discussion of the conditions under which positivity of the splittings is expected. 

In section ^ we use the theory to calculate splittings of a model system and compare to the exactly 
determined values. The agreement is impressive; for the bulk of the splittings, the theory is correct 
to within a few percent. 


2 Extracting the Mean 

We consider double wells with predominantly chaotic dynamics. Our detailed calculations are for 
potentials V(x,y) in two dimensions, though the general results and conclusions apply to higher 
dimensions. A specific example of such a system is provided by the potential, 

V(x,y) = (x 2 -l) 4 + x 2 y 2 + fiy 2 , (1) 

which will be used in section ^ for explicit numerical testing of our results. For positive /r, this 
potential forms two wells, symmetric in x and separated by a saddle at the origin with energy E = 1, 
as illustrated in Fig. [l]. 

In a one-dimensional double well, the energy level splittings A E n vary monotonically from one 
doublet to the next (doublets are indexed by n). In fact, the splitting is a simple, smooth function 
of the mean energy E n of the doublet. In a chaotic multidimensional well, this picture suffers a 
qualitative change — while in any given energy range there is a typical scale for splittings, individual 
splittings vary quasi-randomly from one doublet to the next. The primary purpose of this paper is 
to discuss these fluctuations in the tunnelling rate and to relate them to the properties of individual 
wavefunctions. In order to do this in a coherent way, it is useful first to calculate the typical size of 
splittings in any given energy interval and to use this to rescale the calculation so that fluctuations 
are always of the order of unity. That is the object of this section. We define a simple function of 
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Figure 1: Here are potential contours of a double well with chaotic dynamics. This potential, corresponding 
to Eq. ([!;) with y = 0.1, will be used in numerical investigation of the results. Also represented is the complex 
bounce orbit crossing the barrier in imaginary time —itq and defining the tunnelling action Kq. Linearisation 
of dynamics about this orbit defines a complex monodromy matrix Wq . The classical data Kq and Wq together 
determine the mean splitting function f 0 (E). 

energy which predicts the mean value of the splittings. This is then used to define a rescaling of the 
splittings whose fluctuations are essentially universal. 

In H we introduced the splitting-weighted density of states 

f(E) = Y J ^E n 8(E-E n ) 

n 

Knowledge of this spectral function gives a complete specification of the tunnelling properties of the 
system. A semiclassical expansion of f(E) was found in terms of complex periodic orbits crossing the 
barrier. Here we limit ourselves to the dominant term, 

\ e ~Ko/h 

71 \J (—l) d_1 Ao-det(H / o - 

which furnishes the average properties of the splittings. It is calculated from what is commonly called 
the “instanton” or “bounce” orbit, illustrated in Fig. [lj. This is a complex orbit which crosses the 
barrier with imaginary momentum and real position in an imaginary time t = —iro and with an 
imaginary action S = iK$] it can be understood in terms of a real orbit of the system in which the 
potential is turned upside down. For two-dimensional potentials which are symmetric in y, the orbit 
simply runs along the x-axis bouncing between the two energy contours. Wo is then a complex 2x2 
monodromy matrix linearising motion transverse to the orbit, which can be put in the form 

Wo = ( W0 ~ iv ° ) . ( 4 ) 

V IVo Wo) 




( 2 ) 
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The real numbers (wq, Uq,Vq,Wq) are simply the matrix elements of the real monodromy matrix in the 
inverted problem. Equality of the diagonal elements follows from a symmetry of the orbit under time- 
reversal (this structure is altered somewhat in the presence of a magnetic field and will be discussed 
more fully in section |ib3|). In any case, even though Wq is complex, the amplitude term det (Wo — I) 
is real, as are Wo’s eigenvalues. Finally, d is the dimension of the system and Xa is simply a sign 
fixed by the nature of the symmetry operation a that underlies the double well. When a is an 
orientation reversing transformation of configuration space, Xc = +1 and in the orientation preserving 
case, x. a = — 1- In particular, = +1 in the case where a is a reflection in one coordinate, as in the 
problem used to illustrate the calculation here. In practice, it can be determined simply by demanding 
that the argument of the square root be a positive number. 

The function fo(E) yields the average splittings in the sense that if f{E ) is averaged within a 
sufficiently large window, the contributions of other complex periodic orbits to f(E ) are suppressed. 
This means that the average value of splittings, when measured in units of the mean spacing between 
doublets, is (po(E n )AE n ) = fo(E). We denote by po(E) the average or Thomas-Fermi density of 
states within one well (computed for a given parity in y if appropriate) — its inverse is the mean 
spacing between doublets. Explicit numerical verification that fo(E) describes the mean behaviour 
can be found in [Q]. 

Note that fo(E) typically increases rapidly with energy and therefore so do the typical values of 
the splittings. It is more convenient to work instead with normalised splittings y n , defined by dividing 
by the local average, 

S "~ h(En) ”■ ( ) 

These have average value of unity, by construction, and can usefully be compared between different 
parts of the spectrum, or even between different systems. 


3 Derivation of the Matrix Element for Splittings 


In investigating fluctuations in the tunneling rate, one option is to include, in the formalism alluded 
to in the previous section, complex orbits with real parts to their actions. Analysis of this kind was 
initiated in |i| and a systematic enumeration of the most important such orbits will appear in [ 191. 
This approach is most powerful on large energy scales, where relatively few complex orbits suffice 
to reproduce extremely detailed collective behaviour of tunnelling rates. When the tunnelling rate 
associated with an individual level is required, however, practical implementation of this program can 
be daunting due to the large number of long orbits which must be calculated to reach the energy scale 
of a typical level spacing. For this reason it is of value to pursue an alternative approach, as we do 
in this section, which steers a course between purely semiclassical and fully quantum calculation to 
specify individual tunnelling rates. 

Our discussion follows the methodology of AK ||] and WH Q, whereby each splitting is expressed 
as an overlap between approximate states localised in each of the wells. As discussed in WH, for chaotic 
systems, these local approximations are not found using semiclassical approximations. Instead, it is 
assumed that the wavefunctions in the wells are found through some numerical diagonalisation or 
possibly through some random matrix theory modelling if we are interested in statistical properties 
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of the splittings. Once the wavefunction has been constructed in the wells, the splitting is computed 
by extending the wavefunction under the barrier using semiclassical approximations to the Green’s 
function. 

We give in this section an explicit derivation of the main formula. The purpose of this is two-fold. 
First, it is important for practical implementation to emphasise certain boundary conditions placed on 
the complex classical trajectories used in the formalism, having to do with the path that time takes in 
the complex plane. In addition, we would like to emphasise certain coordinate-invariant properties of 
the calculation which might be useful for application to more general problems. While broadly similar 
to the AK and WH results, the final form of our calculation differs from theirs in certain technical 
aspects, mostly having to do with exercising a freedom as to which parts of the wavefunction are used. 


3.1 An abstract Herring’s formula 


The starting point is Herring’s formula [12]. In our derivation we assume |L) and | R) to be symmetric 
and antisymmetric combinations of odd and even eigenstates |±) which are respectively localised in 
the left and right wells. (While deriving the the splitting for an individual doublet, we will drop the 
doublet label n). We denote the symmetry operation which relates the two wells (typically a reflection 
or an inversion) by the symbol a and call U(a) the corresponding unitary operator. We then have 
| R) = U(a)\L). If we write the coupled Schrodinger equation as 


H\L) = E\L) — (AE/2) |i?) 

H\R) = E\R)-(AE/2)\L), (6) 

then the eigenenergies are = E^fAE/2 and the splitting A E = E~ — E + is positive when the even 
state is at a lower energy than the odd state. Now let P be any Hermitean operator such that P\L) ~ 0 
and P\R) ~ | R). A typical choice would be to let P represent multiplication of the wavefunction by 
the characteristic function ,\s( x ) of the region to the right of some section £ (of codimension 1) which 
separates one well from the other. The precise form will not play an important role here. Applying P 
to each of the equations above and taking matrix elements, we find 


(R\PH\L) = E(R\P\L) - (AE/2)(R\P\R) 

(L\PH\R) = E(L\P\R) — (AE/2) (L\P\L). (7) 


Subtracting the complex conjugate of the first from the second we get 


(L | 


P,H 


I R) = ^ ((R\P\R) - (L\m) . 


( 8 ) 


The calculation is exact up to this point. We now take advantage of the approximations (L\P\L) ~ 0 
and (R\P\R) « 1 to obtain, 


A E 


2 (L\ P,H | R). 


(9) 


The error in this approximation is governed by leakage of the wavefunctions of |L) and | R) across 
£. This effect is exponentially small so the above approximation is very accurate; its error is almost 
certainly dwarfed by that of the next step, which is to substitute approximations for |L) and \R), to 
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be discussed in the next subsection. Also, we will henceforth dispense with using approximation signs 
for the splittings and similar quantities, it being understood that any expression using semiclassical 
arguments is an approximation. 

This explicit expression for A E, which is independent of representation and holds for any form of 
the Hamiltonian, is a generalised version of Herring’s formula. By choosing P suitably, it might even 
be used for calculating dynamical tunnelling splittings, where the division between the localised states 
is in phase space and not in configuration space. The important point is that it relates the splitting 
to an overlap of the localised states that is restricted to the forbidden region of phase space, where 
the Weyl symbol of P is nonconstant. The form used in 

AE = fi2 X ds ’ (io) 

is a special case obtained when H = p 2 /2 + V(q) and P is chosen to represent multiplication by 
the characteristic function xs(x) discussed above [it is a straightforward consequence of the relation 
[x(q)iP 2 ] = ih(x'(Q)P + PX'{o)) which holds for any function %(</)]. Here ds is a surface element along 
£ and the normal derivative d/dn is directed to the side of £ containing \R). In this case we get 
an integral along a lower-dimensional surface because xs( x ) rises from 0 to 1 as a discontinuous step 
function. Had we used instead a smoothly rising function, the overlap would have appeared as a fully 
d-dimensional integral, where d is the dimension of the system, restricted to the region around £ where 
Vxe(x) 0 0. Other possibilities might be possible if P mixes momentum and position operators. 

Matrix elements of the kind shown in Eqs. @ and (|l0[ ) will repeatedly appear during the course of 
the calculation, so it is useful at this point to introduce explicit notation for them and point out those 
aspects which are important. We assume in general that we have a related pair (P, £), where £ is a 
section and P is a Hermitean operator whose Weyl symbol rises from 0 to 1 on passing from one side 
of £ to the other. In Herring’s formula £ is in the forbidden region, but in later incarnations, we will 
have conventional sections in the allowed part of phase space. Defining the anti-Hermitean operator 


As 


P,H 


we construct a sectional matrix element between two states |0} and \ip) defined as, 


( 11 ) 


(0 //e 0) = (01 As |0). 


( 12 ) 


Note that (0 //s V*)* = — (0 //s 0 )- This change in sign might usefully be interpreted as a change in 
orientation of £, or the replacement of P by I — P. The sectional overlap defined in Eq. (12) will 
later in the calculation produce matrix elements in the Poincare section representation introduced in 
[10] to define the transfer operator. It is also similar to notation used in j2j]. When H is of kinetic- 
plus-potential type and P corresponds to a characteristic function xs(x) as above, the matrix element 
becomes the standard surface integral 


//s 0) — 


h 2 


ds 0 


dip 

dn 


d0* 

dn 


0 


(13) 


familiar from Green’s identity, the unit normal being directed to the side of £ where xs is 1. Herring’s 
formula, 

AE = 2{L//^R) (14) 
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now gives the splitting as a sectional overlap of the left and right localised states, albeit in a region of 
phase space where dynamics is complex and the wavefunctions exponentially small. 


3.2 An abstract Green’s identity 


Herring’s formula is formally interesting but not of direct use without further analysis since explicit 
calculation of the wavefunction in the forbidden region can be as difficult as finding the splitting itself. 
Following ||, we simply assume the wavefunction to be known in the allowed region and then use 
Green’s theorem to continue it into the forbidden region, using a semiclassical approximation for the 
Green’s function constructed around complex classical trajectories. 

Let us define sections El and Er = oEl cutting through the left and right wells respectively, as 
illustrated in Fig. |2|. Let V be the region between them. Let Vj, u t e r be a larger region enclosing V and 
G(x, x', E) a Green’s function which satisfies Schrodinger’s equation (with delta-function source) as 
long as x and x ' are within V(, u ter- We will not impose on G(x, x. 1 , E) that it satisfy the Schrodinger 
equation, with boundary conditions, outside of V) m ter■ It may diverge at infinity, for example. In fact, 
it will be important later that it not be the global Green’s function because it is then well-defined 
at quantised values of E. which is where we will want to use it. This is analogous to the use of 


free space Green’s functions in the boundary-integral method in billiards [13]. In our case, we will 
use a semiclassical approximation for G(x, x',E) in which trajectories are discarded after they leave 
Louter ■ The resulting function will then satisfy the necessary equations inside Kuter and has the added 
advantage that, with a judicious choice of V^uten only a finite number of trajectories contribute for 
a given argument. We will further suppose that we only ever want to reconstruct the wavefunction 
within some region Vi nner enclosed by V. Let x(x) be a function which is zero outside I4 ute r and rises 
to unity passing through El or Er and is constant by the time V[ nner is reached. The geometry is 
sketched in Fig. |j. 








Figure 2: A schematic representation of the three domains I4 ute r, V and Vinner- 











That G(x, x',E) is a Green’s function inside V^ uter amounts to the following in abstract notation: 


XG(E - H)x = x(E - H)Gx = x\ 


(15) 


where G and x are the abstract operators corresponding to G(x,x\E) and x( x )- Note that x does 
not have an inverse. Let \ijj) be any state whose wavefunction satisfies the Schrodinger equation in 
14 u ter with energy E. so that x(E — H )|^) = 0. Then we can write the following two relationships, 


xG(E-H)xW) = X 2 \ 
xGx(E-H) |V>) = 0. 


Taking the difference gives, 


XG x,H 


= X 


Dehne the anti-Hermitean operator A = x, H in analogy with Eq. (0. Then, provided x 


ViT, 


V’(x) = (x|GA|^) = (x| G//tp). 


(16) 


(17) 


is m 


(18) 


The overlap here is as in Eq. (|T^) except that the integration is along two sections E^ and E^ instead 
of the single section E. It is an abstract form of Green’s identity capable of relating the wavefunction 
in the centre of the barrier to its values in the neighbourhoods of E^ and E/j. Note that it was not 
fundamentally important in deriving this relation that x represent multiplication by a function of 
position — it is only necessary that Eq. (0) hold and that |?/>) locally obey the Schrodinger equation. 
This observation is helpful in generalising the result to the case of dynamical tunnelling where we 
would want to deduce the behaviour of the state in one part of phase space from its behaviour in 
another. However, we will continue with the language of tunnelling between wells to avoid confusion 
of notation. 

Let us now consider using this method to extend the wavefunction ^(x) from a neighbourhood of 
T,l to the centre of the barrier. We make two approximations, based on the premise that tunnelling 


effects are required only to leading order in A E. First, Eq. fll8|) is employed even though Vt( x ) is 
not, strictly speaking, an eigenfunction. To justify this we note that inclusion of the off-diagonal term 
in Eq. (||) results in the following correction to Eq. ([IT]), 


A P 

XG x,H \L) = x 2 \L) - — xGx\R) 


(19) 


We claim that the additional term on the right hand side can be neglected when evaluating the 
wavefunction near the centre of the barrier. This will be seen self-consistently, first by ignoring 
the term and, following semiclassical approximation of the Green’s function later, noting that the 
imaginary action underlying xGx\B) is similar to that underlying x 2 \L). Therefore, since A E is small 
compared to other energy scales, the second term can be ignored. The second approximation is that, 
since i/jl ( x ) is small in the right well, the contribution to 'ijjL ( x ) in Eq. ( |l8| ) from the neighbourhood 
of E r is ignored. The validity of each of these assumptions can be explicitly verified in the case of 


WKB approximations in one dimension. Thus we can recast Eq. (18) as 


Vt( x ) = ( x l G//z l L) 


( 20 ) 


9 











where // E denotes a sectional overlap confined to E^. With the obvious analog for ^(x), we are 
now ready to invoke Herring’s formula. 

Denoting, 


G tun = G* 


P,H 


g = gV/ s g, 


( 21 ) 


which is a concatenation of advanced and retarded Green’s functions along the section E, we have, 


AE — 2 (L //s^Gtun //•£ r R)- 


( 22 ) 


The minus sign arises from the conjugation of Vt( x ) following calculation from Eq. ( |20| ) (cf. comment 
just below Eq. ©). Having substituted our continued wavefunctions into Herring’s formula, we get 
an expression for A E involving three sectional integrals. There is one integral through the centre from 
Herring’s formula itself, and two more around the surfaces of section coming from the wavefunction 
extensions. In our notation the central integral becomes hidden in the definition of Gtun- This is 
deliberate because, as will be seen in the next section, this integral can be performed analytically 


within semiclassical approximation. The matrix element in Eq. (22) subsequently involves just two 
integrations, along E l and E R . 

Formally, Eq. (|2^) is essentially the same as Eq. (3.8) of However, the formula is implemented 
differently in our case — Auerbach and Kivelson choose sections inside the forbidden region whereas 
we choose them running through the middle of the wells. In addition, the derivation of Eq. (^2| ) could 
equally well be valid for dynamical tunnelling if the operator P is chosen appropriately. The main 
complication of dynamical tunnelling will come later, when the complex trajectories used to construct 
G tU n need to be calculated, in which case the structures involved are less simple than in the case of 
tunnelling between wells. 

In section ^ we will develop an interpretation for Gtun as the transfer operator for a complex 
Poincare mapping. For conventional real Poincare mappings, the transfer operator is a finite-dimensional 
unitary matrix, quantising the classical first-return map |To| . We will define a complex mapping from 
E r to itself for which Gtun acts like a quantisation in the same way (though it is not unitary in the 


complex case). Using the symmetry to identify the left and right wavefunctions, Eq. (22) is then 
formally like an autocorrelation function. 

For the sake of concreteness, we write down the explicit forms taken by these integrals when the 
sectional overlaps are of the form given in Eq. (|P|). First denote V n = d/dn — d/dn, the arrows 
indicating whether the derivative acts on the function to the left or to the right. Then the central 
integral, defining G tU n(x,x', E) = (x|G t un|x'), is 


Gtun(x,x / ,E) = %- f ds" G t (x,x",E) V n n G(x", x', E) 
2 J s 


and the outer integrals are 


A E = -%- f ds [ ds' ^(x) V n G tun (x, x', E) V n , 
2 J Si J 


(23) 


(24) 


Remember that our convention is that the normals to E^ and E# are each directed towards the 


forbidden region. Eq. (24) is the form given in [f)l] (a factor of 2 appears because we use a different 
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normalisation of the Green’s function). There are actually four terms involved, depending on the com- 

4 —> 4—> 

binations of directions in which the derivatives in D n and T> n i act. Semiclassically, their superposition 
conspires to cancel all but trajectories that pass through T,^ and E/j in fixed senses, as we will see. 

3.3 The semiclassical Green’s function 

The standard semiclassical approximation for the Green’s function G(x, x', E) takes the form of a sum 
over trajectories a which go from x' to x with energy E, 


G(x,x',E) = — 


1 


1 


iti (2irih)( d W 2 
For each orbit, S a = f a p ■ dx is the action and 




JSa/Ji 


(25) 


D n = —- det 


xx 


d 2 S a 

dydy' 


(26) 


assuming a single coordinate x along the trajectory and d — 1 coordinates y transverse to it. Detailed 


discussion can be found in [ 141. The usual presentation of this formula has additional complex phases 
determined by the Maslov indices. In our case we assume that they are accounted for by the choice 
of branch for the square root of D a . We do it in this way so that complex trajectories can later be 
included using the same formula, simply by complexifying the classical data entering into it. 

In real WKB the sum is normally taken over all orbits a whose time of evolution is positive. 
This corresponds to giving E a small positive imaginary part and results in the retarded Green’s 
function G(x, x',E). We will also need its hermitean conjugate [see Eq. (^)]. This is the advanced 
Green’s function and can be interpreted as a sum over all orbits whose time of evolution is negative. 


Tunnelling effects are treated by including complex trajectories [15]. This means that the initial and 
final conditions in phase space can be complex, and that time varies over a contour in the complex 
plane, but the energy is still taken to be real. The complex trajectories contributing to the retarded 
and advanced Green’s functions are respectively those obtained by restricting the time contour to the 
positive and negative half-planes. 

The complete Green’s function for the problem requires inclusion of long trajectories, bouncing 
indefinitely in the wells. On the contrary, we will restrict x and x' to the domain fo U t er and include 
only trajectories which take a direct path from x' to x without leaving it. As long as caustic surfaces 
are avoided, which we assume to be the case, this will result in a function which obeys the defining 
equation for the Green’s function while its arguments remain in Io U ter- Thus it fits the description of 


the operator G used in the abstract Green’s identity Eq. (18). The penalty for the restriction on orbits 
is that the semiclassical approximation will be badly behaved at caustics and will ultimately diverge 
exponentially on the forbidden side of them. This will not concern us as long as the caustics of the 
trajectories we use lie outside of fo U t er - We need only a local solution. In particular, the elimination 
of multiple reflections within a well means that G does not develop singularities as the energy passes 
through energy levels. 

Our particular interest is to evaluate the integral in Eq. (|2^). This requires that we evaluate the 
retarded Green’s function connecting the initial position x' near E/j to a point x." under the barrier and 
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Figure 3: A schematic illustration of the trajectories a and (3 respectively determining G(x",x', E) and 
G^(x, x",E) in the integral defining Gtun- The arrows indicate the velocity directions. Note that we expect 
that the trajectories be generically complex, so that the paths also have imaginary parts (not shown). 


subsequently the advanced Green’s function from x" to the final position x near E/,. We will assume, 
following the restrictions of trajectories to Kuter discussed above, that there is a single trajectory for 
each leg of this journey, denoted respectively by a and f3. Consider a. In a one-dimensional potential 
the trajectory starts at x' with velocity pointing towards the barrier and, after an evolution of time 
along the positive real axis, reaches a turning point. A subsequent time evolution along the negative 
imaginary axis allows this orbit to penetrate the forbidden region with an action that has a positive 
imaginary part. For the second part /?, time continues to evolve along the imaginary axis until the 
turning point on the left hand side is reached. Subsequent evolution of time parallel to the negative 
real axis (appropriate to the advanced Green’s function) allows the trajectory to reach x with a final 
velocity pointing once again towards the barrier. Notice that the trajectory begins and ends at points 
which are symmetric in phase space. It is easy to imagine a deformation of this picture into higher 
dimensions, if x, x" and x' are moved away from a symmetry axis for example. The basic structure 
is sketched in Fig. ||. 

To complete the calculation, a saddle point analysis of the integral in Eq. (|23|) is invoked, following 
substitution of these complex orbits into Eq. (p5|). The saddle point condition specifies that the final 
momentum of a matches the initial momentum of /?, so that a concatenated trajectory 7 = (3 o a is 
defined going from x' to x. It turns out that the ensuing saddle point integration gives an expression 
for Gtun(x, x', E) that is precisely of the form given in Eq. (p5|) except that the single orbit 7 is used 
instead of the sum over a. This can be verified by explicit calculation but it is a rather standard 
manipulation and we will not go into it further here. 

Before inserting G tun (x, x 7 , E) into Eq. (24) and writing out explicit results, it is helpful at this 
point to specify more precisely the coordinates x = (x, y) that we will use. We assume such a coordinate 
set to be constructed locally around each of the sections E/, and Er so that the positive x direction 
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points towards the barrier and the area element is ds = d d ~ 1 y. Under normal circumstances, when 
Y>l = all R , we will go further and assume the coordinate set around to be the symmetric image of 
the coordinate set around E#. An important aspect of this convention is that, once a coordinate set 
is fixed on E^, the orientation of the coordinate set on E^ depends on the nature of a. For example, 
in two dimensions the frame on E^ is oriented oppositely according to whether a is a reflection about 
an axis or an inversion through a point. This seemingly innocent difference will have an enormous 
effect on the nature of the final result. When calculated with respect to these frames, the matrix 
linearising motion about 7 is hyperbolic in the case of reflection symmetry and inverse hyperbolic in 
the case of inversion symmetry — we will find as a direct consequence of this that the splittings are 
always positive in the former case whereas in the second they are sometimes negative. This will be 
discussed fully in section ||. Until then, however, our formulas will not depend overtly on whether we 
have inversion or reflection as a symmetry as long as we keep to our conventions for ( x,y ). 

Now define the following modification of Gt un (x, x', E), 

G{y,y',E) = hyjxyxty G tun (x L (y),y, x R (y'),y', E). (27) 

We have scaled out by the velocities which will be reabsorbed into the wavefunctions below; this mirrors 
the development in [l 0 |]. Q will always be evaluated with arguments on E^ and £r, as appropriate, so 
the x-dependences are determined by the implicitly-defined functions xl(jj) and x R (y'), and can be 
suppressed. Taking advantage of the approximations, —ifidGt nn /dx ~ x 7 G tun and —ihdGtun/dx' « 
— x'^Gtun, we may rewrite Eq. ( |22| ) as follows, 


A E = 2 h 


d y [ d y £{y,y') G(y,y',E), 


(28) 


where, 


£{y,y') 


(v x + x 7 ) ^ L (x L {y),y) 

(v'x + x 1 ^ il> R (x R {y'),y') 

2\/x; 

2y / x!) 


(29) 


Here we let v x denote the projection of the velocity operator —iKV onto the direction defined by the 
coordinate x. If we ignore the fact that each of x 7 and x 7 depend on both initial and final position, 
£(y,y') has the appearance of a product of wavefunctions, 


£(y,y')=rfl (y) ^ R {y'). 


(30) 


Each of the wavefunctions here is meant to be identified with a quantity within square brackets in 
the equation above. The abuse of notation inherent in Eq. (^) is justified by the fact that we will 
find that the dominant contribution to the integral comes from a region of the (y, ?/)-plane small on 
classical scales, so that x 1 and xl, can in practice be replaced by constants. 

The arrows on the wavefunctions indicate directions along which flux is predominantly flowing, 
which arises as follows. The surface E# defines two surfaces of section in phase space, E_r and E_r, 
corresponding to the two senses of x. If we look at a phase space representation such as the Wigner 
function, the combination of x 7 and the velocity operator in Eq. (^) leads to enhancement of 
around T, R and suppression around Eh (and equivalently for ipL,)- Therefore we can take the arrows 
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to indicate left- and right-going wavefunctions. The particular combination seen in Eq. ([}(]) arose 
because we chose retarded Green’s functions to extend both of the states into the forbidden region. 
Had we chosen other combinations of advanced and retarded Green’s functions, x 1 and x' would 


have entered Eq. (29) with different signs, leading to different combinations of left and right going 
wavefunctions. It is clearly preferable that this choice be made so that there is symmetry between 
the initial and final coordinates. In particular, once we have chosen between advanced and retarded 
Green’s functions for wavefunction continuation, the same choice is applied to both wavefunctions 
in the matrix element for A E. In phase space, this means that the splitting is calculated from the 
behaviour of states either near £ r and TiR= ct £l, or else near £ r and T,r= ctYiL- 


Im(t) 



Figure 4: The paths taken in the complex time plane by the segments a and (3, and their concatenation 7 . 
Deformation of these contours is allowed as long as branch points are avoided. 

A consequence of this symmetry is that, taking into account the complex conjugation of ^(x) 
in the matrix element, the time evolution of 7 is directed along opposite directions of the real axis 
before and after barrier penetration (Fig. |j). Therefore there is cancellation in the real time along 
the path 7 and the net time of evolution is directed predominantly along the negative imaginary axis. 
For the special case of the periodic bounce orbit 70 , this cancellation is complete and the total time 
—itq is pure imaginary. In computing these trajectories, it is tempting to deform the time contour 
so that it flows directly down the imaginary axis, forgoing the real time segments. In general that 
is dangerous, however, unless one can guarantee that branch points are not crossed as the contour 
is deformed. Besides, the contour defined above gives trajectories whose position coordinates remain 
predominantly real throughout the evolution, which has obvious intuitive advantages. 

To reflect the near-imaginary time evolution, let us use the notation, 

K (y,y', E ) = -iSj(x L (y),y,x R (y'),y',E) (31) 

to denote the “Euclidean” action of the path 7 , suppressing x as usual. As with the velocities x and 
x', the quantity K is predominantly positive but is somewhat complex in general, implying that S is 
dominantly positive imaginary but with some real component. An explicit expression for the modified 
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Green’s function is then. 


G(y,y', E ) = 


(2irh)( d b/ 2 


Xo-det - 


d 2 K 

dydy' 


,-K/n 


(32) 


Note that this has a different ^-dependence than (25) and also has no velocity prefactors. The sign 
Xa is the same as appeared in Eq. (||) and accounts for any reorientation of the local y coordinates 
relative to the global cartesian coordinates with which we started. Except for the fact that it uses 
complex trajectories, the construction of Q(y,y',E) is precisely analogous to that of the Bogomolny 


transfer operator [10]. Likewise, the arrowed wavefunctions can be interpreted as states in the related, 
semiclassically-defined Hilbert space. This interpretation will be developed further in section |j. 

It is at this point that our presentation differs most markedly from previous developments. The 
discussion leading to Eq. (|32| ) is invalid when the end points of 7 are too close to the energy contour 
E = V. In § a careful treatment of this problem is given in which a uniformised semiclassical 
treatment of the Green’s function, using Airy functions, is developed. Auerbach and Kivelson g 
choose to avoid the turning point singularities by taking the surface to be well inside the forbidden 
region. Our approach is similarly to assume and E# to be well separated from the energy contours 
around the dominant trajectories, but on the allowed side. In particular, this means that we must avoid 
the natural temptation to continue the wavefunctions from as close to the barrier as possible, but does 
have the advantage that the ensuing calculation is much simpler (later we will also appeal to symplectic 
invariance to argue that it should be possible to compute in a Poincare section symplectically rotated 
in phase space so that it remains well-defined even at the turning point). It is also advantageous only 
to consider the wavefunction in the allowed region and to put all exponentially small quantities and 
the attendant complex structure into the Green’s function and quantities derived from it. 

As a function of y' and y, the real part of K(y, y', E), which controls the magnitude of Q(y, y', E), 
is minimised when the orbit intersects the Poincare sections with real momenta. This condition is 
fulfilled by the periodic bounce orbit 70 , for which the net time of evolution is imaginary and K is 
real. Recall that this is also the orbit from which /o is determined in periodic orbit theory. It can 
be argued quite generally that the segments of this trajectory corresponding to real time evolution 
take place in real phase space, connecting to a tunnelling segment of imaginary time evolution and 
complex coordinates at a turning point where x = 0 if there is time-reversal symmetry, or Imx = 0 
more generally [fT7i|. 17]. In the potentials treated here, the bounce orbit is on a symmetry axis for 


which y = y' = 0. When its arguments move away from 70 , the Green’s function then decays in 
magnitude over a length scale of 0(Ti 1 ^ 2 ). 

We will assume that the dominant contribution to the splitting matrix element comes from this 
neighbourhood. Approximating the amplitude of the Green’s function by its value for 70 , Eq. ( |2il ) 
becomes 

AE = 2hg(y 0 ,y 0 ,E) [ dyf dy' £(y, y') e~‘ K ^ n (33) 

Aj, 

where IC(y,y') is the quadratic form giving the second variation of Ji(x, x. 1 , E) about the intersection 
of 70 with the sections (for which the coordinates are denoted yo). At the same level of approximation 
we may choose £{y,y') to have the decoupled form (]3^) by taking x and x! to be the values for the 
central orbit 70 . IC(y,y') is calculated from the complex monodromy matrix W defined by 70 , which 
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has the following block form [each block having dimension (d — 1 ) x (d — 1 )], 

W=( W ' '“V (34) 

V IV w ) 

where u and v are hermitean (real if d = 2). That W is of this form follows from symplecticity and 
the property 

W~ x = W*, (35) 

which itself is a consequence of the fact that taking the complex conjugate of the path of 70 in phase 
space gives its time reversed partner, up to the symmetry operation a. The form in Eq. (fT~i|) is less 
restricted than in Eq. (Q) because 70 includes real orbit segments at either end. Standard manipulation 
using generating function conditions obeyed by the action now gives, 


£(jm/0 = \{v y')( a b 



where, 


a = wu 


-1 


and 


b = —u 


-1 


(36) 


(37) 


are respectively symmetric and hermitean. We note that y and y' are in general (d — l)-dimensional 
and the various sub-matrices are also of dimension d — 1. In the case d = 2 they are scalars. Notice 
that, as an operator kernel, JC(y,y') is formally hermitean (it is only formally so because we have 
not yet specified the Hilbert space). This hermiticity is a consequence of the inversion symmetry 
underlying Eq. ( j35[) — i.e. that complex conjugation amounts to time reversal — and is not confined 
to the quadratic truncation. To aid practical implementation, we note that the amplitude in Eq. ( fi^ j 
is det [— d 2 K/dydy'] = detu -1 . 

It is important to note that the integration in Eq. ( |33| ) must be treated with care, since the 
wavefunction will typically vary significantly on the length scale over which the kernel of this matrix 
element decays. In particular, if we were to implement some sort of WKB approximation of the 
wavefunctions (such as a collective approximation using Green’s functions, for example), the center 
of our expansion would not be a saddle point of the whole integrand — further saddle point analysis 
would require that we revert to the full expression for G(y, y ', E ) in that case. In order for the quadratic 
phase expansion to be valid, the integration must be interpreted as an exact sum. Even with this 
limitation, Eq. (|33|) is useful because it can in practice be evaluated following a relatively painless 
numerical diagonalisation to get the wave function, as we will see. 

In any case, we have reduced calculation of the splitting to evaluation of a simple matrix element. 
Either Eq. (28) or its approximation (33) could now be used to determine the splittings semiclassically. 
However, we choose to postpone this step because there remains additional structure which has not 
yet been addressed. In particular, there is an appealing canonical invariance which becomes manifest 
when the theory is reformulated in a phase space representation, as we now discuss. 
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4 A Phase Space Formulation and Scarring 


The matrix element for A E developed in the previous section is conveniently interpreted in the Weyl 
formalism. The Weyl symbol of any operator O is defined by, 

W 0 ( q,p) = J ds e ip s / n <q + s/2|0|q-s/2>. (38) 

and for a state \i/j) the Wigner function W^(q, p) is defined by substituting O = \ip)(ip\/(2Trh) d . The 
expectation value of O can then be calculated as follows, 

{'ippl'i/j) = J dqdp W^(q, p) Wo(q,p). (39) 

We would like to do the same for the matrix element in Eq. (|33|). 

Because we deal with integrations which are restricted to sections (which we will assume to be 
defined by fixing x) it is convenient to define the partial Weyl symbol 

Wg(C, E) = f ds e ipyS / n Q(y + s/2, y - s/2, E) (40) 

and the partial Wigner function 

Wn (0 = 27 Th j ds e ipyS / n £{y + s/2, y - s/2), (41) 

where we denote, 

c =U)- (42) 

We have taken the opportunity here to reintroduce the index n labelling the different doublets. We 
will often refer to ]/V n (£) as the Poincare-Wigner function, which we interpret as a function on the 
oriented section E_r. Canonically invariant expressions will be given below for Wg. 

These functions will usually be evaluated with x and x’ restricted to sections and we therefore sup¬ 
press these variables. These two quantities are then thought of as functions of the section coordinates 

y and p y . In general these are of dimension d— 1, in which case the s integral is of the same dimension 

and the term ip y s in the exponential should be treated as a dot product. Notice that the choice of 
normalisation for yw (£) is not conventional. This will simplify the appearance of the final result. 
It should be remarked that the partial Weyl function could be defined for any operator defined (at 
least semiclassically) on the Poincare section by replacing Q with the operator of interest. Similarly 
we can replace £(y,y') by ^(y)^* (y 1 ) if we are interested in the properties of a semiclassically defined 
wavefunction tp(y) defined on the section. Finally, we note that these transforms are strictly speaking 
properly-defined only when (x, y) are linear, cartesian coordinates, so that the sections E/j and S^ are 
planes. However, since we will always confine ourselves to semiclassical approximation, and usually 
to integrations confined to the immediate neighbourhood of the orbit 7 , it should not be a problem 
in practice to work with curvilinear sections. 

In analogy with Eq. (|39|), the integrations over E# in Eq. (^ 8 |) can be converted to a phase space 
integral on the section using the partial symbols as follows, 

AE " = i k R ( 2^7 Wg{C E ” ) E, ( 0 - («) 
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To evaluate this we need an explicit expression for Wg(C, E). 

Wg(£,E) is structurally similar to the standard Weyl symbol of a time propagator. An explicit 
semiclassical expression for the latter was developed in [1^]. Since the detailed calculation is easily 
transcribed from that case, we will just write down the results. For a given point £ on Er, Wg(C, E) 
is constructed using the “midpoint orbit”. This is the (assumed unique) orbit of type 7 = (3 o a going 
from C, A on Sr to aC B on S l such that the midpoint is £ = (Ca + Cb)/ 2- We associate with it the 
complex action, 

MC,E) = K(x B ,x A ,E) +ipy(yB ~Va)- (44) 


The contribution coming from the y degree of freedom can be interpreted as the symplectic area in 
the y-Py plane of a chord defined by the trajectory and the straight line connecting initial and final 
points. This is invariant under linear canonical transformations in ( y,p y ). Geometrical interpretation 
is more awkward for the full expression, but canonical invariance within a surface of section is certainly 
maintained. We can now write [18], 


Wg(C,E) 


e -A(£,E)/h 

Vx*det [(W AB TIM 


(45) 


where Wab is the complex symplectic matrix linearising motion about the midpoint orbit in ( y,p y ) 
coordinates and Xa is the usual sign compensating for any reorientation of the local coordinates y 
with respect to the global ones. 

The Gaussian approximation employed in Eq. (|33|) corresponds to replacing Wab by its value W 
on 70 and using the following quadratic approximation for A, 

AC E ) - K 0 » C = (c. j^c) (46) 

where T denotes transpose and f2(£, 77) = Jr\ is the symplectic two-form, J being the unit symplectic 
matrix. This result is derived starting with Cb = ^Ca an d employing some generating function con¬ 
ditions satisfied by K(x b ,xa, E) to compute the second derivatives of A. A more complete discussion 
can be found in |l 8 | (beware, however, that J is defined with the opposite sign there). 

For us the feature of most immediate interest is that A is real for real ( y,p y ). This follows from 
the property W~ l = W* discussed earlier. This inversion symmetry implies that the matrix 

fw -I\* _ W-1-I _ W-I 

\W + lJ W~ l + I W + I v ’ 

is imaginary, which in turn, taking account of the factor of i in Eq. ([h]), implies that A(£,E) is real 
to quadratic order. We are therefore free to write 


A(C,e)^k 0 + c t mc 


(48) 


where the matrix 


M = -iJ 


W-I 
W + I 


is real, symmetric and, it turns out, positive definite. 


(49) 
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If we factor out the function fo(E) representing the mean splitting, the Gaussian approximation 
to the symbol Wg(C, E) becomes 


Wg(C,E) « 7r(2vrn) d - 1 /o(^)ffM(C), 


(50) 


where the Gaussian 


9m{ C) 


v^det M -C T Mg/h 
(' nU) d ~ 1 


(51) 


is normalised so that / d£ gM(C) = 1- Note that we have used the fact that W is related to the matrix 
Wo of section | 2 ] by the similarity transformation W = S~ 1 WqS where S is the real matrix linearising 
motion from Sr to the barrier (see below and Fig. |5|). 

The splitting (43) can then be put in the form, 


AE n = f 0 (E n ) l dCsM(C) VVn(C)- (52) 

JXr 

Alternatively, we can write 

Vn = Po{E n ) [_ dC 9m(C) Wn (C) (53) 

for the rescaled splitting defined by y n = po(E n )AE n /fo(E n ) in section |j. 

Eq. (|5^) is the central theoretical result of the paper. Even though the calculation leading to it 
may seem complicated, this final form is very simply interpreted and implemented. Let us summarise 
the ingredients needed to use it. The first step is to find the bounce orbit which crosses the barrier 
in imaginary time and is responsible for the contribution fo(E) to Eq. (p). This reduces to a trivial 
one-dimensional problem in the common case that it lies in a symmetry axis, but is a straightforward 
numerical task even in the most general case. Quite generally this bounce orbit has a turning point 
P where its coordinates in phase space are real. Take P and its symmetric image crP to be the 
initial and final points of the imaginary-time trajectory, denoted 7 r. Starting at P, one can equally 
consider evolution in real time, defining a real trajectory 7 r going in negative time to £ r . Its image 
7 l = < 77 r goes to aP from £ r. Then the concatenation 70 = 7 r 0 7 r 0 7 r, where the star denotes 
time reversal, is an orbit going from Sr to Si in a net time which is imaginary. The monodromy 
matrix W = S^WoS of 70 defines M in Eq. ©• The final ingredient is the computation of the 
Poincare-Wigner function >V n (£). There are no simple semiclassical theories in the case that the 
dynamics is chaotic. It can be calculated numerically however, with the advantage that the precision 
need not be excessively high. One can also model the wave function from random matrix ensembles, 
and Eq. ( |53| ) will form the basis for quantitative statistical analysis of the splittings in that case. 

In cases of additional symmetry, one encounters the situation where the real extension to 7 r is 
a periodic orbit. In this situation it is common to interpret an exceptionally large accumulation of 
probability amplitude in its neighbourhood as a manifestation of scarring j|, [l 8 |]. Our expression for 
the y n is in such cases simultaneously a “scar-ometer” — it gives a precise measure of the degree to 
which the state is scarred by the periodic orbit. In fact, it is extremely close to being the value of a 
Husimi distribution evaluated on the orbit. It would be precisely so if 3 m(C) were the Wigner function 
of a coherent state. It is not — to be a Wigner function of a pure state requires detM = 1, which 
is not true in our case — but it is very close. In particular, a proper Husimi distribution would take 
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Figure 5: A schematic representation of the structures needed to calculate the splitting. First the imaginary- 
time trajectory connecting the surfaces of section on each side of the barrier must be computed. This defines a 
complex monodromy matrix of the form W = S' _ 1 WoS' where S is the real matrix linearising motion around a 
real trajectory segment. W defines a Gaussian 3 m(C)i approximately supported in an elliptical domain of area 
0(K) in each Poincare section, as indicated by the shaded regions. Overlap with the Poincare-Wigner function 
of the state in this region gives the splitting. 


the same form as Eq. (53) except for a rescaling of the matrix M. The other difference between our 
result and a standard Husimi functon is that our overlap is confined to a lower-dimensional Poincare 
section. The fact that there is a strong correlation between scarring and the tunnelling rate in such 
cases is not at all surprising. Scarred states have a large accumulation of probability along the optimal 
tunnelling route defined by 70 and would be expected to tunnel at a higher rate. The fact that this 
observation can be given a precise quantification is nontrivial, however. 

As a scar-ometer alone, the expression has the following interesting, and potentially useful property: 
our calculation indicates, and explicit numerical calulation in section |7] will verify, that the overlap 
is an invariant of the orbit. That is, it does not matter for which section £/j (he. which value of 
x) we evaluate the overlap, the result is always the same within semiclassical error. It gives a single, 
unambiguous number with which to quantify the degree of scarring of the state. An important element 
in this invariance is that the elliptic region defined in the Poincare section by level curves of C, T M£ 
deforms with the flow as we move from one point of the orbit to another — as a consequence of 
W = S _1 WoiS we have M = S t MqS, where S linearises the flow along the real orbit. In particular, 
as we move further from the turning point at which tunnelling begins, this elliptic domain becomes 
stretched along the stable manifold of the periodic orbit. Clearly, if the evolution is taken too far, errors 
of linearised propagation will grow and the invariance will fail. However for moderate deformations 
invariance is expected and observed. 

Finally, we would like to verify that the result is correctly normalised. To this end we note that 
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the general property of standard Wigner functions, 


'£W^ P)S (E-E n) \ = W*Az3 


(54) 


can be argued (appendix |I|) to have the following analog in the case of Poincare-Wigner functions; 


£W„(0 5(E-E n )} =xe(C). 


(55) 


Here Xe(C) is the characteristic function of the energetically-allowed region in the £-plane [i.e. Xe{C) 
is unity whenever £ is within the energetically permitted domain and is zero otherwise.) The average 
here is meant to be taken over an energy range small on classical scales but large enough to smooth 
out quantum fluctuations. This tells us that the Poincare-Wigner function is expected on average to 
be spread uniformly over the surface of section in the canonical measure dydp y , which is unsurprising 
(note however that renormalistion of the wave function by the velocity as in Eq.(|29|) is a necessary 
ingredient for this to work). A fuller discussion of this, along with certain qualifications, is offered in 
appendix jB[ As a consequence, 

(Y,y n 6 (E - E n )j = po(E). (56) 

which follows from taking the average inside the £ integral of Eq. 
words, ( y n ) = 1 as required. 

We note finally that fluctuations of J2 n Wn (C) $(E — E n ) about the mean defined above can 
be treated semiclassically as a sum over real midpoint orbits returning to Er (l^]. We can combine 
these sums with our matrix element for y n to derive semiclassical expressions for the fluctuations of 
— E n ). In doing so we simply reproduce the expression for f(E) as a sum over complex 
orbits which was found using the trace formula |J. Though less direct, this derivation has technical 
advantages because the long orbits can be effectively treated within real dynamics and ambiguities due 
to the possible intervention of Stokes phenomena etc, normally difficult to control in higher dimensions, 
disappear. 


53) and taking gu (C) out. In other 


5 The Case of Resonances 


Now we outline how the discussion of the previous sections may be adapted to the calculation of 
resonant lifetimes in metastable systems. Many of the formal expressions in this case look identical. 
There are important differences in detail, however, particularly in the boundary conditions placed on 
the imaginary-period bounce orbit which mediates tunnelling. 

To be concrete, we will initially assume a potential with a fairly restricted topography, to be 
relaxed later. As illustrated in Fig. ||, we assume a single, classically bound potential well, with a 
barrier leading to a single escape channel. In quantum mechanics, we find approximate states localised 
within the well. These are associated with poles of the scattering matrix at complex energies 


E n = E„ - 


iTr 


(57) 
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Figure 6: Potential contours of a typical metastable system. The section E c han is defined so that, restricting 
integration to the left of it, the resonant wavefunction is normalised. The resonance width is related to the flux 
across it. The wavefunction is assumed known on the section £ we ii cutting through the well and a semiclassical 
Green’s function is used to propagate out to £ c han- The orbit shown is appropriate to the retarded Green’s 
function and corresponds to type a. 


where T n is the exponentially small resonance width which is the object of discussion in this section. 
(Note that there will also typically be resonances with very large widths associated with resonance 
wavefunctions which have negligible amplitude inside the well. Their widths are not mediated by 
tunnelling so such states do not concern us here.) 

Letting \ijj n ) be the metastable state associated with a particular resonance, we formally write the 
Schrodinger equation, 

H\^ n ) = E^ n ) - (58) 

Being a metastable state, | ip n ) is unnormalised, but it is not difficult to define an operator P so 
that P\i/j n ) is both normalised and an approximate eigenstate. Typically, we contract some region 
A enclosing the well and excluding all but finite portions of any escape channels and let P represent 
multiplication of the wavefunction by its characteristic function. For example, in Fig. || we might 
let A be the region to the left of S c han- We will further assume \i/j n ) to be normalised so that 
(r/> n |P|?/> n ) ~ 1. Applying P to the equation above, taking matrix elements and then subtracting the 
complex conjugate, we get the following modification of the abstract Herring formula, 


Fn — 


P,H |V>n) = 


(59) 


Physically, this relates the resonance width to the flux leaving the well. 

As before, we use a semiclassical Green’s function to extend the wavefunction from the interior of 
the well across the barrier and into the escape channel, allowing evaluation of Eq. (|5^) along £ c han- 
We assume a section £ we u can be defined in the well such that a single complex trajectory connects it 
with E c h an while remaining free of caustics in the region enclosed between them (or at least this should 
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be so where the dominant contributions to integration arise). In other words, this single trajectory 
defines semiclassically a Green’s function G which satisfies Schrodinger’s equation in the region where 
Green’s theorem demands it. Such a trajectory will be labelled by a and is illustrated in Fig. |b[ 
Notice that a completely crosses the barrier and has real-time segments before and after the crossing. 
Thus it looks superficially like the orbit 7 of the previous section. The important difference is that, 
in the typical case where we use the retarded Green’s function, the real time evolution of a is along 
the postive axis before and after the crossing. 

These real parts will be cancelled when we integrate the Green’s function against its hermitean 


conjugate, following evaluation of the diagonal matrix element in Eq. (59). To see this more explicitly, 
we note that calculation of the resonance width using the sequence of the previous section leads us to 
define the following tunnelling operator, 


4u„=iGt //^J, 

in terms of which the resonance width is, 


(60) 


— (V’n //s^n^tun //s^nV’n)- ( 61 ) 

The sectional overlap restricts integration to the neighbourhood of E we ii and is defined in obvious 
analogy with the previous sections. Except for the factor of i, Gt U n looks formally identical to its 
splitting counterpart defined in Eq. (^). As in that case, the real time evolutions coming from G and 
G' will cancel — following saddle point integration, Cr tun is approximated with a single orbit 7 = /3 o a 
obtained by concatenating an orbit a of positive real time evolution with an orbit j3 of negative real 
time evolution. 

Despite its formal similarity with that of the previous sections, it is important to stress that the 
single orbit 7 from which the present Gtun is constructed has very different boundary conditions 
imposed upon it. In particular, it crosses the barrier twice , beginning and ending on E we u. It will 
be obtained by deforming an imaginary time orbit which is periodic, not pseudoperiodic as in the 
previous section (in any case, there is no longer a symmetry with respect to which pseudoperiodicity 
could be defined). 

Following detailed pursuit of this program, we arrive at the following Poincare matrix element for 

r„, 

r n = h f d y f d y' G(y,y',E) £ n (y,y'), (62) 

J ^well J ^well 

where £ n (y,y') is defined analogously to Eq. (f29|). The Poincare Green’s function Q(y,y',E) takes 
exactly the same form as given in Eq. (|32|), except that the imaginary action K (y, y'. E ) is defined by 
the double-bounce orbit appropriate to resonances as described above. Note that this is in spite of 
the factor i in Eq. (^), which is compensated by different phases which arise on passing through the 
barrier on the way to E c h an , as can be verified in standard one-dimensional WKB. 

Having performed the saddle point manipulation leading to Q(y,y',E), the real time segments 
connecting the outer turning point of imaginary evolution to E c h a n can be dispensed with, and in fact 
Dchan disappears from the discussion, just as the intermediate section did in the previous sections. 
Beforehand, however, the detailed saddle point analysis is more straightforward if S c h an is chosen to 
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^well 


Figure 7: For resonances, the orbit 7 takes the path shown. As usual, the arrows show the velocity direction 
during the real time segments. 


the right of the barrier. If it is taken under the barrier as in the previous section, the requirements 


vplaced on P leading to Eq. (59) are still met, but the dominant contribution to the integral along 
E c han is awkwardly hidden behind direct-path contributions which are formally greater but which 
cancel in a detailed analysis. The best approach is simply to avoid the whole issue by choosing E c h an 
as shown. 

As with splittings, the natural next step is to expand K(y,y',E) quadratically about the max¬ 
imum of its real component. This condition is given by a periodic orbit of imaginary period (with 
real extensions). There are obvious analogs to the matrices W and M defining the quadratic phase 
expansions in the double well case, and these share the symmetry properties derived from W* = IT -1 . 
Finally, we can transform to a phase space representation using the partial Wigner functions as in 
section [|. 

To present the final result, we note that the resonance-width-weighted density of states 


f(E) = J2^nS(E-E n ). 

n 

is approximated in the mean by the term 


(63) 


ME) = 


e -K 0 /h 


2vr 


\J (—l) rf_ 1 det(ll / o — I) 


(64) 


defined by the bounce orbit, in analogy with Eq. (|64|). Then the rescaled resonance width 

Po{E n ) „ 


Un — 


MEn 


(65) 


has average value unity and is given by a formula formally identical to Eq. (53). 









Finally, we note that it should be obvious how to modify the discussion when different topographies 
arise. For example, in the case of multiple channels one would simply sum over the different bounce 
orbits and let the optimal one dominate. (The same applies to the case of symmetric double wells 
when there are multiple tunnelling routes connecting the wells.) 


6 General Structure of the Result 

This section consists of two parts. In the first a formal interpretation using transfer operators is 
developed for the result, and in the second positivity of the splittings is related to hyperbolicity of the 
symplectic matrix W. 


6.1 A formal interpretation using transfer operators 

The matrix element we have derived for splittings and resonance widths can be formulated very 


concisely in terms of the Bogomolny transfer operator and the structures surrounding it 11C ]. This 
is defined in terms of the real dynamics within a given well, so let us restrict ourselves to that case 
initially. For convenience of presentation we will use the notation developed for the case of splittings 
in symmetric wells — subsequent adaptation to resonance widths in metastable wells will be obvious. 
Therefore we start with the oriented Poincare section E?? and denote the usual first return map by 
F : £i?—>£??. One supposes now that there is a Hilbert space TL(Tir) quantising Si? and then a 
unitary transfer operator T : 7f(Si?) —► TL( Si?) is constructed within semiclassical approximation 
whose classical limit is the symplectic map F : Si?—►Si?. Matrix elements of T are given by a 
formula like Eq. |3^) except that the real orbits and actions defined by F are used. Wavefunctions in 
W(Si?) are related to proper wavefunctions by a two-step process: first the part of the wavefunction 
corresponding to traversal of the section in the correct sense is extracted, and then the amplitude is 
modified by the square root of the transverse velocity. This corresponds neatly to the construction of 
the arrowed wavefunctions in Eq (p9|). In particular, we will interpret a directed wavefunction (ip r (v) 
for splittings and the obvious analog for resonance widths) as the wavefunction of a particular element 
| n) of TL(Tir)- Finally, in the transfer operator formalism proper eigenfunctions are related to the 
existence of eigenfunctions of T with unit eigenvalue — thus we expect T|n) = |n). Note that a precise 
quantum-mechanical definition of these objects has not yet been given in the completely general case. 
This need not concern us, however, because we are using them for the purposes of interpretation only 
and ultimately have well-defined formulas such as Eq. (^) to fall back on. 

So far the discussion has been without reference to tunnelling. We assume these structures have 
been calculated in standard semiclassical approximation within one well, ignoring exponentially small 
effects. These are now incorporated using the tunnelling Green’s function construction of the preceding 
sections. The Green’s function Q(y,y',E) will be interpreted as the semiclassical quantisation of a 
complex mapping F : E/?—>Ej? which is constructed by using orbits of type 7 to map Ej? into itself. 
To understand the distinction between this map and the usual real first-return map F : E??—>Ei?, it 
is useful to think of a one-dimensional double well. Up to symmetry operations, this has two periods 
in the complex time-plane; a real one corresponding to oscillation within a well and an imaginary one 
under which any initial condition is taken to its symmetric partner on the other side of the barrier 
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(with complex intermediate coordinates even if the initial and final conditions are real). The maps 
F and F can then be interpreted as deformations of these periodic motions into multidimensional 
dynamics. It is interesting to compare their behaviours under conjugation, 


F* = F 

F* = F~\ ( 66 ) 


which is what one would expect of real and imaginary time evolutions respectively. Note, however, that 
F and F are not equal-time mappings. The bounce orbit 70 defines a real fixed point of F. However, 
F is in general a completely complex map — other initial conditions, even if real, are mapped to 
points with complex coordinates. Thus when we write F : we must interpret Sfl as being 

a complex manifold. In principle we should expect to be able to extend the definition of F for initial 
conditions arbitrarily far from 70 by analytic continuation on Y,r , though there will presumably be a 
rather complicated branch-point structure in general. In practice semiclassical approximation requires 
that we consider F only in the immediate neighbourhood of 70 , which we assume to be free of branch 
points and singularities. 

Just as there is an operator T quantising F, we can construct an operator T quantising F — its 
matrix elements are given by the modified Green’s function G(y, y',E) in Eq. (|32|). It is not unitary 


like T but Hermitean. In analogy with Eq. ( 66 ) we compare their properties of inversion as follows, 


T x = T 

T x = T- 1 . (67) 


where for any operator we denote with a f the inverse of its Hermitean conjugate (to see why we choose 
this as the quantum analog of complex conjugation, consider the operator exp[— izH\ for complex z 
and Hermitean H and define complex conjugation to be that operation which sends z —► z*). 

We can now formally write the rescaled splitting as, 


Vn = 


(n | T | n) 

Tr T ' 


( 68 ) 


noting that the Weyl symbol of T is given in Eq. (45) or, in quadratic phase approximation, by 


Eq. (51). It should be noted that the state | n) is not normalised conventionally. Formally, it can be 
defined to be the state for which the Wigner function, defined conventionally according to Eq. (filj) 
with d replaced by d — 1 and coordinates restricted to the surface of section degrees of freedom, is, 


Wn(C) = P 0 (En) Wn (0- (69) 

In practice, calculation of \n) boils down to Eq. (^), or some version semiclassically equivalent to it, 
along with appropriate redefinition of the normalisation. Eq. (^) indicates that | n) is normalised so 
that the average value of this Wigner function is unity or, more generally, such that its coefficients in 
a generic basis for H(E_r) average to one. 

While Eq. (| 6 ^) can be used as it stands, it is considerably easier to implement if a Gaussian 
approximation is used, as in Eqs. (|33|) and (^). This amounts to substituting for T the quantisation 
of the linear symplectic map W instead of that of the full nonlinear map F. Before writing this 
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explicitly, it is useful to be more specific about the structure of W. Define R to be the symplectic 
transformation which inverts an eigenvector of W if the corresponding eigenvalue is negative and leaves 
it invariant otherwise. Then R commutes with W and RW is hyperbolic, i.e., has only positive real 
eigenvalues. We then implicitly define a real positive definite symmetric matrix B through 

W = Re~ iJB , (70) 

That B defined in this way is symmetric follows from the fact that W is symplectic and that it is real 
follows from Eq. (|35|). Positivity comes from more specific dynamical properties of the system and is 
discussed further in the next subsection. In other words, we consider W to be generated by evolution 
under the positive quadratic Hamiltonian, 

K C) = \c t bc (7i) 

for the imaginary time At = —i, followed by a reorientation R of certain eigenvectors to allow for 
negative eigenvalues. 

The Gaussian approximation now corresponds simply to the substitution, 

T' = 77exp[— h/h]. (72) 

for T. In the exponent, h is a positive-definite hermitean operator which is quadratric in the operators 

a. ~ T * 

y and p y and whose Weyl symbol is h(£). We will write h = £ BC,/ 2 where £ is a vector of operators 
defined in obvious analogy with Eq. ([42|). U is a unitary matrix quantising R — its sign is fixed by 
demanding that its Weyl symbol be positive. Up to normalisation, 5 m(C) is the Weyl symbol of this 
operator, as can be checked explicitly. 

We note finally that in two degrees of freedom the spectrum of T 7 is the geometrically decaying 
sequence 

Spectrum [ T 1 ] = {(±l) k e^ k+1 /2)ha/h } ^ = {|A|“§, (73) 

of powers of the eigenvalue A = ±e“ of W which is larger in magnitude, a = [detl?] 1//2 being the 
frequency associated with real-time evolution under h(Q. In higher dimensions the spectrum of T’ 
consists of products of such powers of eigenvalues. 

6.2 Reflections, inversions and negative splittings 

In certain situations, one expects that the y n ’s should systematically be positive. This is certainly the 
case for resonance widths, for example. There is also often an expectation in the case of symmetric 
double wells that the even state should be at a lower energy than the odd state, as is the case in 
one-dimensional wells and in the system illustrated in section It is not always so, however. For 
example, when the symmetry underlying the splitting is a complete spatial inversion a : x —> —x, 
negative splittings are routinely encountered for excited states, as discussed in appendix ^]. The 
following question now arises. When is the matrix element underlying y n inheherently positive and 
how can we relate this positivity to to structural and dynamical properties of the problem? 
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The answer is straightforwardly encoded in Eq. (|f^) and can be summarised by the following 
statement, 

Assertion: In Gaussian approximation, T is a positive definite operator whenever W is hyperbolic, 
that is, whenever W has only positive eigenvalues. In particular, this is expected to be the case for 
resonance widths and for splittings in symmetric double wells when the symmetry is reflection in one 
cartesian coordinate. When W has negative eigenvalues, so does T', and we expect both positive and 
negative values for the y n ’s. This is the case, for example, in a symmetric double well whose symmetry 
is a complete inversion. 

Hyperbolicity of W is to be expected for a complex orbit crossing a barrier in a potential and computed 
with respect to a rigidly translated basis (see below). When nonhyperbolicity arises, it is because the 
coordinate frame with respect to which the final displacement is measured has undergone a reorienta¬ 
tion with respect to the rigidly translated basis. In the case of symmetric double wells, such a frame 
transformation would be induced by the application of a to map the coordinate system of one section 
onto its symmetric partner as in the discussion preceeding Eq. ([27|). In the case of metastable wells, 
there is no such transformation and W is always hyperbolic. 

To understand why we expect W to be hyperbolic in the absence of frame rotation, consider the 
special case of a potential barrier with a symmetry axis, so that 70 is quasi-one-dimensional, as in 
Fig. 0. Then W is the solution of the differential equation, 

^ = -iJ[H"] ± W, (74) 

dr 

where [H"] j_ is the Hessian matrix of the Hamiltonian in cordinates transverse to 70 . We expect the 
potential, and therefore [H"\ 7 , to be positive definite underneath the barrier and it then seems clear 
that the solution of Eq. (|74| ) would be of the form e~ lJB , with positive real symmetric B , as in the 
previous subsection. We claim without proof that this is the case and that it also holds in the more 
general case when 70 is not on a symmetry axis. 

It should be acknowledged that, while positivity of T' is certainly a necessary condition for the 
positivity of the y^s it cannot yet be taken as sufficient. It will be seen in explicit implementation in 
the next section that it is often the case that there is significant cancellation in the matrix element, 
leaving a very small result, and corrections could in principle tip the balance and result in negative 
splittings. We do not believe this to be the case in practice, but have not rigourously proved otherwise. 

It is interesting to note how positivity arises in the Wigner-Weyl formalism of section [|— irrespec¬ 
tive of the hyperbolicity of W, the splitting is given by the overlap of a positive Gaussian 5 m(C) with 
the Poincare-Wigner function, and at first sight it is not obvious why and when this overlap would be 
negative. In broad terms, the solution is as follows. If <7m(C) is sufficiently narrow (or equivalently 
if M is sufficiently large) then one can see that it might be supported in a region where the Wigner 
function is negative and a negative overlap might emerge. On the other hand, if <?m(C) is very broad 
(so M is very small), we might expect the overlap to wash out any negative regions and only positive 
overlaps to emerge. Unsurprisingly, the smallness or largeness of M is linked to the hyperbolicity of 
W. 

To simplify the discussion, let us restrict ourselves to d = 2. We can then give a very precise 
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criterion for M in order that the overlap be positive. This is summarised by the following statement, 


Assertion: The overlap of qm (C) with an arbitrary Wigner function is always positive precisely 
when det M < 1, in which case <7 m(C) the Weyl symbol of a positive definite operator of the form 
exp[— C, B£/2h\. In the limiting case det M = 1, 5 m(C) is the Weyl symbol of a pure state and, when 
detM > 1, it is the Weyl symbol of an operator of the form U exp[— Q B£/2h] where U is a unitary 
operator quantising the inversion R = —I. 


The first part of this statement comes simply from calculating the Weyl symbol of exp[—£ B£/2h], 
which is proportional to, 


exp 


tanh(a/ 2 ) 


C t bc 

ha 


(75) 


where a 2 = det-B, and equating its exponent with that of <7m(0- A solution can be found precisely 
when detM < 1 (note that det [tanh(a/2)B/a] = tanh 2 (a/2) < 1). The marginal case detM = 1 
is not directly relevent to tunnelling calculations, but was included for completeness — it would 
correspond to the limit A —> oo. The case detM > 1 is easily seen to correspond to the insertion of 
an inverse hyperbolic matrix W in Eq. |4^. Note that replacing W with — W in that equation leads to 
the following replacement for M, 


W —► W' = -W =>• M —> M' = 


(76) 


In particular, det M = 1/det M, so det M > 1 corresponds to an inverse hyperbolic matrix as claimed. 

~T - 

It can also be verified by explicit calculation that if <?m(C) is the Weyl symbol of exp[—£ B£/2h\, 
then gM'{ C) is the Weyl symbol of U exp[— £ B£/2h\. In other words, <?m'(C) is the Moyal product of 
9m(C) with the Weyl symbol of U. 

Finally, we note that in the case d = 2, there are two possible symmetries for a double well. 
Reflection through an axis defines a hyperbolic W and positive splittings, whereas inversion through 
a point defines an inverse hyperbolic W and produces a mixture of positive and negative splittings 
(appendix [(]). 


7 An Explicit Implementation 

To test the results, we will work with the potential defined in Eq. ([!]) of section [2], and illustrated in 
Fig. Classical motion is constrained to one or other of the wells for energies in the range 0 < E < 1 
and we expect nearly degenerate doublets in the quantum spectrum. This is a convenient choice for 
the numerical calculation of quantum spectra since the potential is polynomial in x and y and the 
Hamiltonian is then banded in a basis of harmonic oscillator eigenstates, allowing large bases and 
therefore accurate diagonalisation. In addition, the classical dynamics is largely chaotic over a wide 
range of energies. 

In a previous publication j9j], we discussed tunnelling rates in this potential for the case fi = 0. 
This choice has the disadvantage that the potential does not have a generic behaviour at the top of 
the barrier; instead of a saddle point, one finds a degenerate ridge seperating the two wells. This leads 
to nongeneric behaviour of tunnelling rates for states just below this critical energy. A corrective 
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analysis is certainly possible but would be specific to that case and not particularly interesting in its 
own right. We therefore prefer to avoid this situation by choosing p > 0. On the other hand, large 
values of /1 induce large regular regions in phase space and so we limit ourselves to values of p that 
are relatively small, specifically p = 1 / 10 . 

The classical motion is then predominantly chaotic as long as the energy is not too small (typically, 
E > 0.3). A typical Poincare plot is shown in Fig. || for E = 0.7. The section is defined by plotting 
y and p y whenever x = 1 and p x < 0. The real periodic orbit which connects with the bounce orbit 
intersects this Poincare section at the origin. Although some small regular islands are present, they 
are hardly visible and the dynamics is dominantly chaotic. It is an important feature because we can 
then expect the wavefunctions to be ergodically spread throughout the space so that they can “access” 
the central tunnelling route 7 discussed above. For nonchaotic systems, wavefunctions are typically 
strongly localised in phase space and the previous analysis would be relevant only to a subset of the 
states, assuming the real extension to the bounce orbit to be unstable. If it is stable, the region around 
it supports eigenstates whose energies and splittings can be extracted semiclassically [9|] without any 
reference to the wavefunction at all. The presence of chaos will also be important to justify various 
arguments necessary in a statistical analysis of the splittings @]. 



Figure 8 : Classical dynamics is illustrated in the form of the x = 1, p x < 0 Poincare section for E = 0.7 and 
p = 0.1, which is the case used for numerical calculation. The section is almost completely chaotic, with no 
visible islands. 

In order both to invoke this formalism and to test it, we need to diagonalise the quantum problem 
numerically. On the quantum-mechanical level, there is an additional parameter which appears, 
namely Planck’s constant. The classical mechanics is independent of Ti while it does depend on the 
choice of energy. For this reason, the analysis is much cleaner if we hold the energy fixed at E and find 
those values of Ti for which it is an eigenvalue — i.e. perform a numerical quantisation of fi. This can 
be done by considering Schrodinger’s equation to be a generalised eigenvalue equation for fi. Defining 
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( 77 ) 


the momentum squared operator p 2 /2 = h 2 T, we have 

{E -V)\^ n ) = nlf\^ n ) 

This is readily solved in a two-dimensional harmonic oscillator basis using standard numerical routines 
which make use of the fact that T is a positive-definite operator. In what follows we will analyse the 
spectrum of the reciprocal quantity 

1 (78) 


Q = 


h 


The even and odd spectra are very close but there are small tunnelling splittings A q n between the 
respective states, just as for the energies. 

The only remaining issue is that the theory developed was for splittings in energy so we would 
like to convert our (/-splittings into energy splittings. This can be done using first order perturbation 
theory. Each doublet consists of one state which is even with respect to reflections in x and another 
which is odd; the even state always being observed to have the lower energy. The same property 
holds as a function of q\ the even states have slightly smaller values than the odd states. This can be 
understood if we think of each state corresponding to a parametric curve in E — q space. The doublets 
are then pairs of very closely spaced curves in this space. Because T is a positive-definite operator, a 
first order perturbation calculation shows that dE/dq < 0. From this it follows that if the even state 
of each doublet has a smaller value of energy at fixed q , then it also has a smaller value of q at fixed 
energy. First order perturbation theory then gives 


A E n = - T (V’n|7 1 |V ; n)Ag rl 

Qn 


(79) 


The matrix element is readily calculated, since the wavefunction i/j n is known and so, given the q 
splitting, we can readily determine the energy splitting. Any approximation for i\) n that is more 
accurate than semiclassical error is acceptable. In practice this would mean performing a numerical 
calculation with a basis confined to one well. The relation (|79| ) is valid only to leading order in A q n , 
however this level of approximation is consistent with earlier arguments. 

We have tested the formalism in describing the splittings for the choice E = 0.7, /x = 0.1 and even 
y-parity. We have also worked predominantly with the surface of section defined by p x < 0 and x = 1, 
though some results will also be presented for the choice x = 1/2 (these sections are respectively 
referred to as E a and S^). For the former choice we have that the real time of evolution for the orbit 
7o is to = 0.93 while the imaginary time (i.e. for the bounce orbit) is — itq = —1.237 Therefore, 
we start at E# = E a with y = p y = 0, x = 1 and p x given by energy conservation (with p x < 0). 
We integrate these initial conditions for a time to, which brings us to the barrier. Following this we 
integrate for a time — itq, which carries us across to the left-hand well. Finally, we integrate for a time 
—to, bringing us to E^, with the final velocity pointing once again towards the barrier. The resulting 
point is at y = p y = 0 and x = — l and p x > 0, which is the mirror image in phase space of the initial 
point. While doing these integrations, we also keep track of the complex monodromy matrix W. We 
then construct the matrix M as given by Eq. ( |49| ) and thereby construct the Gaussian function on the 
Poincare surface defined in Eq. (ID- 


As a final step, we evaluate the partial Wigner function Eq. (41) of the wavefunction i^ n . This 
can be done efficiently by expressing the state in a position representation, using the expression (p9|) 
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Figure 9: Some explicit examples of Poincare-Wigner functions, obtained by quantising h for the parameter 
values illustrated in the classical surface of section of Fig. || (negative values are represented by white space). 
Cases (a), (b) and (f) have exceptionally large splittings and correspond to states scarred by the periodic orbit 
at the center. Their Wigner functions are concentrated along its stable and unstable manifolds, shown as the 
heavy curves emerging from the origin. State (c) has an exceptionally small splitting, and (d) and (e) are 
“typical”. Also shown in each case is the ellipse along which the Gaussian 5 m(C) has fallen by a factor e -1 / 2 
its area is nh/2 [detAfp ss 5 h. This elliptic structure is defined by the complex monodromy matrix W and 
is distinct from the hyperbolic structure of the invariant manifolds associated with the real matrix S which 
linearises dynamics about the center. This structural dichotomy emerges even though the eigenvalues of both 
matrices are on the real axis. 


to determine the semiclassical wavefunction defined on the section and then performing a fast fourier 
transform to get the phase space representation. Since the Weyl transform of the Green function 
above was calculated using the central Gaussian approximation, we do the analogous thing for the 
wavefunction by taking the velocity factors appearing in ( |4d| ) to be independent of y and y' and simply 
equal to the values for the central bounce orbit. A few cases are shown in Fig. || where it is shown 
that states with relatively large splittings are strongly scarred while the other states are not. In fact, 
for the most strongly scarred states, one can even observe amplitude extending out along the stable 
and unstable manifolds of the orbit. This supports the statement that the quantity y n and the phase 
space integral leading to it represent a quantification of the extent of scarring. 

We then integrate this function with the Gaussian function as prescribed in Eq. ({Td|) . One final 
complication is that the symmetry in y necessitates some additional analysis. In particular, states 
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even in y typically have larger splittings than states odd in y simply because they have more amplitude 
on the tunnelling path. This property is handled automatically in the first equation of Eq. (53) since 
the factor >V n (£) will typically be larger for the even states than for the odd. Therefore in evaluating 
the raw splittings A E n , we use the same prefactor /o for both symmetry classes; the one described 
in section |2| However, we will want to normalise by different amounts for the two classes to get y n 
values which separately average to unity and so we need to use parity-specific forms for /q and po- 
Specifically, we use the relation y n = p^(E n )AE n /f^(E n ) to get the normalised splittings, where the 
± refers to the y parity and the means for determining these functions is given in appendix |A]. In 
this way, we finish with normalised splittings which average to unity for both symmetry classes even 
though the average splittings before normalisation are different. 

Again, we stress that this only makes use of the wavefunction in the classically allowed region and 
is not a particularly demanding numerical task. 



Figure 10: Some explicit numerical results, calculated in the Poincare section = £ a defined by x = 1. 
The horizontal axis is just an index labelling the states and the vertical axis shows the rescaled splittings on 
a logarithmic scale. The exact quantum results are shown as circles. There are two separate semiclassical 
calculations, represented by triangles and crosses. The triangles represent the simplest calculation using the 
Gaussian approximation — these are joined by lines to the quantum results so as to guide the eye. Agreement is 
generally good, except for very small splittings, which tend to be overestimated by the approximation, though 
the situation improves somewhat with increasing quantum number. These smaller splittings are often better 
estimated by the second semiclassical approximation, represented by the crosses, in which the exact action and 
amplitude of the complex orbit is used. This is the noncentral analysis referred to in the text. 

In Fig. 0 we show the results of this procedure (triangles), for the first 100 states even in y, 
and make a comparison to results of a completely quantum-mechanical diagonalisation (circles). The 
agreement is rather impressive, especially for the large splittings. For the bulk of the states, those with 
rescaled splittings of order unity or greater, the agreement is typically within a few percent. It should 
be remarked that the agreement does improve systematically as one climbs higher in the spectrum. 

The states with small splittings are anti-scarred — having suppressed amplitude along the central 
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periodic orbit. This is problematic since our Gaussian approximations assume that it is the central 
region which dominates the integral. While certainly true for the bulk of the states, it may not be 
true for the handful of anti-scarred ones. To test this, we went beyond the Gaussian approximation. 
Specifically, we numerically found the complex orbit connecting each pair of points (x 1 , y') and (. x , y) 
in the two wells. Each such orbit has an action, a stability matrix and initial and final x velocity. 
We fed this information directly into Eqs. (29) and (fl2|). From those, we could then determine the 
corresponding Weyl and Wigner functions using Eqs. ( |40|) and (41). These are still real functions 
due to the hermiticity properties satisfied by £ and Q. We then used Eq. (|43|) to determine the 
splittings. The results of this analysis are represented by the crosses in Fig. E The agreement is 
dramatically improved for those states with small y values — as expected. It is also generally improved 
throughout the spectrum, except for the very first few states. It should be stressed that these are still 
only approximate since there are still semiclassical approximations being made in the analysis. Note 
that, this procedure is far more difficult than just using the central approximation since thousands 
of complex orbits must be found rather than just the one trivial central orbit. A compromise which 
seemed to improve the results significantly while still requiring relatively little work was to expand x 
and x' as quadratic forms in y and y' about the central orbit and using these expressions in Eq.(|fJ). 
However, for the sake of brevity we do not include these numbers. It should be mentioned that in the 
semiclassical limit Ti —> 0, we expect that only a vanishing fraction of the states will be so anti-scarred 
as to require this non-central analysis. 

In Table ||, we give a more detailed comparison of results for some selected states — the first 
five and a second sequence of five higher in the spectrum. A comparison is made between quantum 
calculation and of our matrix element evaluated in two distinct sections, E a and Eb as well as with 
using the non-central analysis discussed in the preceding paragraph (for the section E a ). Note that the 
results from E a and E& are different, but consistent within the level of approximation achieved while 
those from the non-central analysis improve significantly the small y values. On the whole, results 
from E a or E& seem more or less equivalent. However, we have noticed that the smallest splittings 
tend to be better reproduced by the E^-calculation. 

A possible explanation is that semiclassical errors increase with increasing real time propagation 
at either end of the barrier crossing. This suggests that it might be advantageous to perform our 
Gaussian overlap of the wavefunction as close as possible to the turning point. As we have discussed, 
the derivation we have outlined is problematic in this limit. However, we expect that the final result, 
being formally canonically invariant, should be equally valid as long as a consistent generalisation of the 
Poincare-Wigner function is employed. If calculation is then done in a Poincare section that remains 
well-defined at the turning point, we expect the main formula to remain valid. The main ingredient 
in such a procedure is to perform a symplectic rotation of the quantities entering main result — since 
the Wigner-Weyl calculus is covariant under symplectic rotation this is certainly possible, though we 
have not implemented it in detail. 

We have verified that the formalism works for other choices of energy, parameter values and y 
parity, although, for sake of space, we do not present the corresponding results. 

Finally, we stress that while a very precise quantum mechanical calculation has been possible for the 
potential in Eq. this in large part due to its polynomial nature and is not at all a straightforward 
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n 

q n 

Vn [QM] 

Vn Pa] 

Vn P&] 

Vn [NC] 

1 

2.361 

0.5475 

0.4628 

0.3801 

0.4408 

2 

5.067 

1.2132 

1.1521 

1.1238 

1.1276 

3 

5.924 

0.2267 

0.3407 

0.2315 

0.1437 

4 

7.896 

1.8310 

1.7773 

1.8864 

1.6533 

5 

8.442 

0.5282 

0.5122 

0.5940 

0.4227 

81 

38.795 

0.0601 

0.0891 

0.0742 

0.0653 

82 

39.235 

0.0050 

0.0132 

0.0055 

0.0048 

83 

39.260 

1.4533 

1.4842 

1.6126 

1.4004 

84 

39.506 

0.5517 

0.6219 

0.5964 

0.5484 

85 

39.839 

5.0215 

5.0546 

5.1767 

4.9092 


Table 1: A detailed comparison of rescaled splittings for the first five doublets, and for another sequence of 
five higher in the spectrum. The third column gives results of an exact quantum mechanical diagonalisation. 
Semiclassical approximations were constructed using Poincare sections at x = 1 (£ a ) and x = 1/2 (£&), results 
of which are respectively shown in the fourth and fifth columns. The final column shows the results using £ a 
but with the non-central analysis. Clearly there is no formal reason to expect agreement for the lowest states, 
and the comparison of the first five should be regarded as a case of experimentation rather than verification. 


task in general. In other potentials, for which equivalently large bases are not feasible, it may well 
be impossible to calculate energy levels to the precision necessary to determine level splittings when 
these are very small. On the contrary, implementation of the matrix element formula will hardly be 
different in terms of difficulty. It is for such systems that the formula will be of greatest help as a 
numerical method. 

Although we have not yet tested the formalism for the case of resonances in a metastable well, this 
is something we intend to do in the near future. 

8 Conclusion 

We have shown that an energy-level splitting or a resonance width can be effectively calculated as 
a matrix element of an explicitly-constructed tunnelling operator, defined on the quantisation of a 
Poincare section and having a Gaussian kernel. This tunnelling operator is completely determined by 
the action of an imaginary-time “bounce orbit” and a linearisation of dynamics around it. Besides 
providing a relatively painless numerical algorithm for calculating splittings (otherwise a delicate task 
numerically), we hope that the formal simplicity of the result will be useful in theoretical analysis, 
especially statistical ||20| . 

Even though the bounce orbit crossing the barrier is complex, it defines a completely real orbit 
in each well, sharing a turning point with each of them in the case of time-reversal symmetry. The 
matrix element measures the weight of the wavefunction about these real extensions — the complex 
nature of the crossing enters only in the definition of the Gaussian kernel, which is constructed using 
the complex monodromy matrix W of the bounce orbit. Though complex, W defines a real positive- 
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definite quadratic “Hamiltonian” h(y,p y ) defined on a Poincare section, associated with which is an 
elliptic structure centered on the real orbit. The behaviour of a Poincare-section Wigner function, 
within an area of 0(Ti d ~ 1 ) defined by these elliptic domains, suffices to determine the splittings (except 
for the strongly anti-scarred ones for which a non-central analysis seems necessary). In particular, in 
the common case where the real extension is periodic, we have a simple interpretation of the tunnelling 
rate as a scarring weight. We stress that computation of the classical data entering in this formula 
is of the same degree of simplicity for all systems or parameter values, while a completely quantum 
calculation of the splitting can easily become impractical if the splitting is extremely small or efficient 
diagonalisation methods are unavailable. 
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A Symmetry Decomposition 

The example we choose to study has an additional reflection symmetry in the y direction as well as 
the reflection symmetry in x which is responsible for the presence of a double well. To evaluate the 
quantities y n we need to normalise by appropriately symmetrised functions p^(E) and f±(E n ) as 
discussed in section [7]. In this appendix we explain how this is done. 


A.l Thomas-Fermi calculation 


We will first work with the integrated densities of states which can then be differentiated with respect 
to either energy or q to obtain the corresponding densities. In fact, we need the densities corresponding 
to a given value of the y parity, TT y = ±1. The problem of symmetry decomposition of the smooth 
density of states has been discussed in various papers both for potential systems |2l| and billiards [221. 
Our case is rather simple since it is a simple parity group and the result is, 


N£(E,q) = ±(N I (E,q)±N R (E,q)). 


(80) 


where Nj(E,q ) and N R (E,q ) correspond to the group elements of identity and reflection through the 
x axis. (The reflection operator through the y axis (i.e. n x = ±1) plays no role in this calculation 
since the axis lies outside of the classically allowed domain.) 

To leading order in q, 

N r(E) = ^JdzQ(E- H(z)) (81) 

where z collectively represents the phase space coordinates (x,y,p x ,p y ). The integral is greatly facil¬ 
itated by remarking that the potential (||) is quadratic in y. Energy conservation then gives 

2 2 

E-U(x) = ^- + ^- + (Xx 2 + p)y 2 , (82) 
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where U(x ) = (x 2 — l) 4 . This is just the equation for an ellipsoid with semiaxes ( a,a,b ) where 
a = \J2(E — U(x)) and b = a/\J 2(Ax 2 + /r); its volume is then 4ira 2 b/3. Therefore in (|8l| ) we have 
only to numerically evaluate the x integral, 


-ZVj(-E) 


- r 

3vr J x _ y Ax 2 + /r 


(83) 


where x± = yliEh 4 . The reflection operator contributes an amount |]2l|| 

= £ / dxd Pa; 0(£ - (p 2 /2 + U(x))) 

= ^£ + dxy/2(E-U(x)). (84) 

To determine the densities of states in E or in q. we differentiate ( |80| ) with respect to the appropriate 
variable. This simple exercise is left to the reader, resulting in quadrature which are readily evaluated 
numerically. It is clear that this only applies for energies below the barrier. Also, we have only 
considered the right well, obviously there is an identical contribution from the left well. However, 
strictly speaking, we are usually interested in the density of doublets and not the density of states. 
There are half as many doublets as states which compensates for the factor of two from not including 
the left well. Finally, these expressions are asymptotic series in fi 2 = 1/q 2 which we have truncated at 
the first term. The first correction to Nj(E) is a constant which is typically small for potentials and 
furthermore has no effect on the density of states at finite energy and q. 


A.2 Stability factor calculation 


The function fo(E) defined in Eq. comes from considering just the pure bounce orbit. In the pres¬ 
ence of a reflection symmetry in y, this orbit lies on the symmetry axis and its amplitudes decomposes 
between the even and odd states in a specific way |23|]. We start by defining the larger eigenvalue of 
Wo appearing in Eq. (;3[) by A. Then we can write the stability factor (with d = 2) as 

1 _ A 1 / 2 

V-det(W 0 - I) ~ A- 1 


A 3/2 A l/2 

a^T + a^T' 


(85) 


The first equation comes trivially from the determinant while the second comes from expanding out 
the factor of 1/(A — 1) as a geometric series in 1/A, separately collecting the even and odd powers 
and then resumming the two series. In any event, it is trivial to see that the second line equals the 
first regardless of how it was derived. 

We now follow reference [23] and identify the first term as giving the contribution to the even 
states and the second term as giving the contribution to the odd states. The exponential factor in 
Eq. (I) remains the same and so we identify 


fo(E) 


1 A 3 / 2 

7T A 2 — 1 


e -K 0 /h 


1 A 1 / 2 

7T A 2 — 1 


e~ Ko/h . 


( 86 ) 
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We observe that the even splittings are, on average, larger by a factor of A. 


B Thomas-Fermi for Poincare-Wigner functions 

In this appendix we outline how to derive the microcanonical background of the Poincare-Wigner 


functions [Eq. (|5^)] from that of the standard Wigner functions [Eq. (54)], which we assume to be 
known 0 - 

The main ingredient is the following transformation rule between ordinary Wigner functions and 
partial Wigner functions, 


W^(x,x',y,p y ) = (2irh) d J dp x e Wx{ - X x ' )/h W$ {{x + x')/2,y,p x ,p y ) , 


(87) 


straightforwardly deduced from their definitions in Eqs. (|39| ) and (40). In particular, if we restrict x 
and x' to a section E defined by x = xq and transform Eq. (FI), we get, 


E W n (x 0 ,xo,y,p y ) S(E - E n )\ = J dp x 6(H(x 0 , y,p x ,p y ) - E) 


= E 


dH 


dp x 


-l 


= Em« 


( 88 ) 


where the sum is over all phase space points on the section E, H = E with coordinates (y,p y ). Note 
that this includes trajectories crossing E in both senses. 

This is not yet the result we need, because we have used the raw wavefunctions instead of the 


renormalised versions in Eq. (3C) to define the partial Wigner functions. The transformations defined 
in Eq. (p^) are satisfactory for the purposes of computing splittings because they conform intuitively 


to the restriction and rescaling inherent in the Bogomolny transfer operator [lOj in the immediate 
neighbourhood of the trajectory 70 . They came naturally from Green’s identity, but we do not claim 
that they represent a clean, global definition of Poincare wavefunction (which should be made without 
reference to a Green’s function for example). Since we are not going to offer such a definition, it is 
not worthwhile entering into a technically detailed discussion of the effects of our particular rescaling, 
and instead we will argue intuitively. 

If we were to expand the wavefunctions in a basis of states with a well-defined semiclassical inter¬ 
pretation (as WKB expressions using Lagrangian manifolds for example), the Bogomolny rescaling is 
semiclassically well-defined — project out the components of the basis states corresponding to crossing 
of E in a particular sense, and rescale the amplitudes by a square root of the transverse velocity. The 
effect on the partial Wigner functions would be to restrict them to oriented Poincare sections and 
then to multiply by the transverse velocity. As a result, we adjust Eq. (| 88 | ) by restricting the sum 
to traversals of E in a fixed sense and multiplying each contribution by v x . We assume that ( y,p y ) 
provide a single-valued parametrisation of the oriented Poincare section. Then, for each ( y,p y ), the 
microcanonical estimate is 1 if ( y,p y ) is in the energetically allowed region of the plane and is zero 
otherwise. This is the content of Eq. (|55|). 

To summarise, we have argued that a coherent, global definition of Poincare-Wigner functions 
should have a microcanical background given by Eq. (F3) . We state without proof that our particular 
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definition, as implied by Eq. (29), satisfies this criterion near the intersection of 70 with El (which 
is the only place it is needed to calculate tunnelling). Global properties, such as the nature of the 
transition from 1 to 0 at the edge of the allowed region, are usefully defined only following a more 
invariant definition, which is not offered. 


C Inversion Symmetry 

As discussed in section ||, positivity of splittings in a double well depends on the kind of symmetry 
underlying them. In two dimensions, for example, reflection-symmetric potentials V(—x,y) = V(x,y) 
are responsible for positive splittings whereas inversion-symmetric potentials V(—x,— y) = V(x,y) 
produce both positive and negative splittings. It turns out that we can see this dichotomy very simply 
in the example considered in Section |7|, which has both symmetries at once. It is instructive, as we do 
in this appendix, to investigate how the formalism depends on which symmetry we use. 

In Section [7] that problem was treated as if the symmetry was reflection in x. As predicted in theory, 
all splittings were then found to be positive. Had we instead defined a to be inversion, however, we 
would have found negative splittings as a result of the following simple observation. A state that 
is called “even” with respect to reflection symmetry can be odd with respect to inversion, and vice 
versa for an “odd” state. In such cases, a splitting will change sign and become negative if we replace 
reflection by inversion as the relevent symmetry. 

To see in detail how this emerges, consider the table |2|. This shows the four classes of states, 
defined by how the transform under the various symmetries. For reflection, we define the splittings to 
be the differences between the energies of class 2 and those of class 1 or between the energies of class 
4 and those of class 3, since in each case this is the odd energy minus the even energy where odd/even 
is defined relative to the a x operation. If instead we define the splittings relative to cr x a y , then the 
table tells us that we should continue to use the differences between classes 2 and 1 but should now 
take the differences between the energies of class 3 and those of class 4. In other words, what we were 
previously referring to as the odd y parity doublets are now responsible for negative splittings. 



I 


Gy 


1 

+ 

+ 

+ 

+ 

2 

+ 

- 

+ 

- 

3 

+ 

+ 

- 

- 

4 

+ 

- 

- 

+ 


Table 2: The four classes of state and how they transform under the various symmetries. 

In the calculation described in the text, this is evident if we look, for example, at Eq. (p4|), but 
having defined iPl(x, y) = iPr(—x, —y) rather than 4>l(x, v) = ^r(—x, y). Since the states have definite 
y parity, this implies that there is simply an overall minus sign for the odd parity states and no such 
sign for the even parity states. How does this sign manifest itself in the final results, such as Eq. ( j53j ) 
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or Eq. (p8|)? The answer is essentially that when we switch from a x to o x cr y as the symmetry, we 
should change the sign of W, to account for a reorientation of axes, and the positive-definite Gaussian 
9m(C) is replaced by a nondefinite one, fjM'iC)- Recall that the frame in which W is calculated is 
defined on El to be the symmetric image of that on E r, and has opposite orientation according to 
whether we define a = a x or a = o x a y . 

An added complication, which has already been alluded to, is that the the matrix Wo defining the 
function fo(E) similarly changes sign. This makes perfect sense. Since some of the splittings are now 
taken to be negative, we expect the average splitting value to be less than before. Specifically, the 


fact that the coordinate y has changed sign implies that Eq. (||) becomes [24 ] 


HE) = - 


o-Ko/n 


7T v/det(W 0 + 7) 


= L e ~Ko/h 

IT 


A 3/2 A l/2 


A 2 - 1 A 2 - 1 


(89) 

(90) 


where Wo is the monodromy matrix when the symmetry is a x and A its larger eigenvalue. The second 
line follows from expanding out the determinant assuming d = 2 along the lines discussed around 


Eq. (85). The two terms in brackets have the interpretation of being the even and odd y parity 
contributions but now being subtracted rather than added, just as we expect. 

Note that this change in fo does not affect the magnitude of the raw splittings obtained from 
Eq. (|5^) . The modified value of fo is compensated by the modified matrix M so that the splittings are 
the same with the exception that some of them are negative. However, the normalised splittings are 
changed. This is clear. If they should average to unity but some of them are now negative, then they 
must all be rescaled to make this work out. This is implicit in Eq. (^) where the modified matrix M 
entering the factor gM is not compensated by the prefactor as in (|5^). 

Throughout this discussion, we have been comparing to the the picture in which we used x reflection 
to define the splittings. However, we may not always have this choice. For example we could add 
a small term to the potential such as exy which preserves the inversion symmetry but breaks the 
reflection symmetry. If e is small enough then we do not expect the splittings to change very much. 
Negative ones will mostly stay negative. So the possibility of negative splittings is generic and not 
some pathology unique to this problem. Another more interesting interaction which would maintain 
the inversion symmetry while breaking the reflection symmetry would be to add a uniform magnetic 
field. It is really that possibility which motivates this appendix. 

Also, this discussion has been specific to two dimensions. One could imagine in higher dimensions 
that the symmetry involves reflecting some of the coordinates on the section El but not all of them. 
Then the discussion at the beginning of the appendix still applies if we change the sign of those 
variables which are reflected and keep the sign of those which are not. It is too tedious to work 
through all possibilities in a unified manner and we refrain from doing so here. 
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